Программирование и научные вычисления на языке Python/§18
Используя математические правила из теории вероятностей можно вычислить вероятность того, что произойдет какое-то событие, например, как это случается в задачах, вероятность достать из шляпы шар какого-то цвета. К несчастью, уже при небольшом изменении условий теоретический расчет вероятностей может зайти в тупик и приходится использовать численные методы. К таким методам относится и метод Монте-Карло. Принципы моделирования методом Монте-Карло Представим, мы провели N'' экспериментов и результат каждого эксперимента случайное число. В этих ''N экспериментах некоторое событие случается M'' раз. Оценка вероятности тогда ''M/N и она становится более точной при стремлении числа экспериментов N'' к бесконечности (заметьте, что дробь при этом не обращается в ноль, поскольку и число ''M также стремится к бесконечности). Математический метод, заключающийся в использовании большого числа генерируемых случайных чисел, получил название метод Монте-Карло. Этот метод оказался чрезвычайно полезным в науке и промышленности, там где случайное поведение нельзя не учитывать, или, например, там где и не подразумевается наличие случайных чисел — в случаях сложного интегрирования. Бросание костей Пусть у нас спрашивают: какова вероятность при броске костей, чтобы по крайне мере на двух из них были шестерки? Итак, эксперимент в том, что кидаются четыре кости, нужное событие - не менее двух кубиков с шестью точками. Собственно, это нам и нужно записать еще раз перефразировав на Python: import random as random_number N = int(raw_input('Enter number of experiments: ')) M = 0 for i in xrange(N): six = 0 r1 = random_number.randint(1, 6) if r1 6: six += 1 r2 = random_number.randint(1, 6) if r2 6: six += 1 r3 = random_number.randint(1, 6) if r3 6: six += 1 r4 = random_number.randint(1, 6) if r4 6: six += 1 if six >= 2: M += 1 p = float(M)/N print 'probability:', p Более общая задача — ведь мы можем задать и не только число экспериментов, а, например, число бросаемых костей (ndice) и минимальное число одинаковых кубиков (nsix). Таким образом, мы получаем более общий и даже более короткий код, а более общие решения гораздо проще модифицировать далее: import random as random_number import sys N = int(raw_input('Number of experiments: ')) ndice = int(raw_input('Number of dice: ')) nsix = int(raw_input('Number of dice with six eyes: ')) M = 0 for i in range(N): six =0 for j in range(ndice): r = random_number.randint(1, 6) if r 6: six += 1 if six >= nsix: M += 1 p = float(M)/N print 'Probability:', p Для малых вероятностей число M'' мало и аппроксимация ''M/N не очень хорошая аппроксимация, требуется большое число экспериментов, которые для стандартного модуля обходятся дорого, то есть долго. Конечно же, мы знаем, что при нехватке скорости надо перейти к более быстрым векторным операциям. Векторизация состоит в том, что мы создаем двухмерный массив случайных чисел, в котором число строк определяется числом экспериментов (бросков), а число столбцов — числом испытаний (костей): eyes = random.random_integers(1, 6, (N, ndice)) Следующий шаг — посчитать число нужных нам событий в каждом эксперименте. Для того, чтобы программа работала быстрее мы должны постараться избежать циклов. В предыдущем уроке мы использовали прием создания массива по выполнению условия, этот же прием мы используем ниже. Далее в полученном массиве суммируем элементы строк, а это нули и единицы, и в случае если число единиц равно или больше требуемому числу шестерок, то считаем, что событие произошло: from numpy import random, sum N = int(raw_input('Number of experiments: ')) ndice = int(raw_input('Number of dice: ')) nsix = int(raw_input('Number of dice with six eyes: ')) eyes = random.random_integers(1, 6, (N, ndice)) compare = eyes 6 nthrows_with_6 = sum(compare, axis=1) # суммирование по столбцам - элементам строки (axis = 1) nsuccesses = nthrows_with_6 >= nsix M = sum(nsuccesses) p = float(M)/N print 'probability:', p И подсчет при этом для большого числа экспериментов, происходит существенно быстрее. Шары из шляпы В шляпе 12 шаров: четыре черных, четыре красных и четыре синих. И мы хотим написать программу, которая достает из шляпы три случайных шара. Шары у нас отличаются цветом, при этом эти цвета неизменны, значит их удобно представить в виде кортежа, а саму шляпу в качестве списка строк: colors = 'black', 'red', 'blue' # (кортеж строк) hat = [] for color in colors: for i in range(4): hat.append(color) Чтобы достать шар: import random as random_number color = random_number.choice(hat) print color Но нам нужно достать несколько шаров. При этом мы не возвращаем их обратно в шляпу, а это значит что элементы (шары) удаляются из списка hat. Три способа: 1) использовать hat.remove(color), то есть достали, например, синий шар - значит одним синим шаром в списке меньше, 2) выбирать шар по случайному индексу, а затем по нему же его удалять через del hatindex, 3) то же самое, но более коротким способом через hat.pop(index): def draw_ball(hat): color = random_number.choice(hat) hat.remove(color) return color, hat # или: def draw_ball(hat): index = random_number.randint(0, len(hat)-1) color = hatindex del hatindex return color, hat # или: def draw_ball(hat): index = random_number.randint(0, len(hat)-1) color = hat.pop(index) return color, hat balls = [] for i in range(n): # n - число доставаемых шаров color, hat = draw_ball(hat) balls.append(color) print 'Got the balls', balls Продлим наше решение на более конкретную задачу, чем просто доставание из шляпы шаров: какова вероятность достать два и более черных шара из этой шляпы. В таком случае мы можем проделать N'' экспериментов и подсчитать сколько ''M раз мы достали не менее двух черных шаров и найти вероятность такого события как M/N. Один эксперимент для нас состоит в создании нового списка hat (перемешивании шаров), доставании шаров и подсчете числа черных. Последнее может быть легко подсчитано с помощью метода списков count, то есть в нашем случае hat.count('black'). Тогда оставим одну понравившуюся нам функцию draw_ball(hat) и допишем нашу программу следующим образом: import random as random_number def draw_ball(hat): color = random_number.choice(hat) hat.remove(color) return color, hat def new_hat(): colors = 'black', 'red', 'blue' hat = [] for color in colors: for i in range(4): hat.append(color) return hat n = int(raw_input('How many balls are to be drawn? ')) N = int(raw_input('How many experiments? ')) M = 0 for e in range(N): hat = new_hat() balls = [] for i in range(n): color, hat = draw_ball(hat) balls.append(color) if balls.count('black') >= 2: M += 1 print 'Probability:', float(M)/N Ограничение роста населения В Китае уже в течение нескольких лет супружеским парам разрешено иметь только одного ребенка. Однако успех в проведении этой политики ограничивается рядом факторов. Одна из проблем в том, что семьи чаще оставляют сыновей, которые могут помочь семье, и отсюда все сильнее растет доля мужского населения. Отсюда альтернативная политика в том, что семьям разрешается завести дочь пока не родится сын. Мы можем промоделировать эту ситуацию и посмотреть как будет расти численность в таком обществе. Поскольку мы работаем с большой популяцией, сразу зададимся тем, что поиск решения будем искать в векторизованном варианте. Представим, у нас есть n индивидуумов, назовем их parents, которым мы придаем некоторые равномерно распределенные случайные значения. Далее из статистических данных мы знаем долю мужчин (male_portion), и всех индивидуумов, у которых их случайное число меньше male_portion, считаем мужчинами и присваиваем значение 1, а у кого выше — женщинами и присваиваем значение 2. Чтобы код было удобнее читать, введем константы MALE=1 и FEMALE=2. Наша задача в том, чтобы посмотреть как меняется массив parents для каждого нового поколения при учете двух указанных политик. Итак, создаем исходный массив: from numpy import random, zeros r = random.random(n) parents = zeros(n, int) # массив нулей размером n MALE = 1; FEMALE = 2 parents< male_portion = MALE parents>= male_portion = FEMALE Число потенциально возможных супружеских пар это минимум от числа мужчин и числа женщин. Однако, даже если мы посчитаем, что все эти пары сложатся, то все равно только часть из них дадут детей. При учете политики «одна семья — один ребенок» получаем: males = len(parentsparents MALE) # число мужчин females = len(parents) - males # число женщин couples = min(males, females) # число потенциальных супружеских пар n = int(fertility*couples) # пары, родившие ребенка, fertility - рождаемость # следующее поколение r = random.random(n) children = zeros(n, int) children< male_portion = MALE children>= male_portion = FEMALE Код, написанный для рождении нового поколения будет нужен при каждой новой генерации. Поэтому эти инструкции более естественно заключить в соответствующую функцию: def get_children(n, male_portion, fertility): n = int(fertility*n) r = random.random(n) children = zeros(n, int) children< male_portion = MALE children>= male_portion = FEMALE return children При условии более мягкой политики, назовем ее «политикой одного сына», родители могут рожать детей, пока у них не появится сын. В рамках нашей оценочной задачи будем упрощенно считать, что они так и делают: # сначала обычное children = get_children(couples, male_portion, fertility) # в случае, если дочь: daughters = children FEMALE while len(daughters) > 0: new_children = get_children(len(daughters), male_portion, fertility) children = concatenate((children, new_children)) daughters = new_children FEMALE Накопленный опыт мы можем представить в виде общей функции анализа популяции для обеих политик: def advance_generation(parents, policy='one child', male_portion=0.5, fertility=1.0, law_breakers=0, wanted_children=4): males = len(parentsparents MALE) females = len(parents) - males couples = min(males, females) if policy 'one child': # у каждой пары только один ребенок: children = get_children(couples, male_portion, fertility) max_children = 1 elif policy 'one son': # каждая пара продолжает рожать детей до тех пор, пока у них не будет сына children = get_children(couples, male_portion, fertility) max_children = 1 daughters = children FEMALE while len(daughters) > 0: new_children = get_children(len(daughters), male_portion, fertility) children = concatenate((children, new_children)) daughters = new_children FEMALE max_children += 1 # кроме того часть граждан нарушает закон и имеет столько детей, сколько хочет (wanted_children) illegals = get_children(int(len(children)*law_breakers)*wanted_children, male_portion, fertility=1.0) children = concatenate((children, illegals)) return children, max_children Законченная программа с заданием исходных данных тогда будет выглядеть так: from numpy import random, concatenate, zeros MALE = 1; FEMALE = 2 def get_children(n, male_portion, fertility): n = int(fertility*n) r = random.random(n) children = zeros(n, int) children< male_portion = MALE children>= male_portion = FEMALE return children def advance_generation(parents, policy='one child', male_portion=0.5, fertility=1.0, law_breakers=0, wanted_children=4): males = len(parentsparents MALE) females = len(parents) - males couples = min(males, females) if policy 'one child': # у каждой пары только один ребенок: children = get_children(couples, male_portion, fertility) max_children = 1 elif policy 'one son': # каждая пара продолжает рожать детей до тех пор, пока у них не будет сына children = get_children(couples, male_portion, fertility) max_children = 1 daughters = children FEMALE while len(daughters) > 0: new_children = get_children(len(daughters), male_portion, fertility) children = concatenate((children, new_children)) daughters = new_children FEMALE max_children += 1 # кроме того часть граждан нарушает закон и имеет столько детей, сколько хочет (wanted_children) illegals = get_children(int(len(children)*law_breakers)*wanted_children, male_portion, fertility=1.0) children = concatenate((children, illegals)) return children, max_children N = 1000000 male_portion = 0.51 fertility = 0.92 law_breakers = 0.06 wanted_children = 6 generations = 10 # начнем с "идеального поколения" родителей: start_parents = get_children(N, male_portion=0.5, fertility=1.0) parents = start_parents.copy() print 'one child policy, start: %d' % len(parents) for i in range(generations): parents, mc = advance_generation(parents, 'one child', male_portion, fertility, law_breakers, wanted_children) print '%3d: %d' % (i+1, len(parents)) parents = start_parents.copy() print 'one son policy, start: %d' % len(parents) for i in range(generations): parents, mc = advance_generation(parents, 'one son', male_portion, fertility, law_breakers, wanted_children) print '%3d: %d (max children in a family: %d)' % (i+1, len(parents), mc) Интегрирование Одним из самых ранних применений генерируемых случайных чисел было вычисление интегралов. Пусть мы генерируем равномерно распределенные между a'' и ''b случайные числа x1, ..., xn, тогда аппроксимация решения будет найдена как Этот метод обычно и называется интегрирование по Монте-Карло. Мы можем представить выражение (18.1) в виде небольшой программы: import random as random_number def MCint(f, a, b, n): s = 0 for i in range(n): x = random_number.uniform(a, b) s += f(x) I = (float(b-a)/n)*s return I Обычно, достаточная точность метода задается больщим числом n'', поэтому векторизованная версия будет более удобной: from numpy import * def MCint_vec(f, a, b, n): x = random.uniform(a, b, n) s = sum(f(x)) I = (float(b-a)/n)*s return I Рассмотрим интегрирование методом Монте-Карло на примере простой линейной функции ~f(x) = 2x + 3 , пределы интегрирования — от 1 до 2. Было бы интересно посмотреть как метод справляется с решением задачи для различных ''n. Оценку произведем следующим немного измененным MCint методом: def MCint2(f, a, b, n): s = 0 I = zeros(n) for k in range(1, n+1): x = random_number.uniform(a, b) s += f(x) Ik-1 = (float(b-a)/k)*s return I Заметим, что k' изменяется от 1 до n'', в то время как индексы ''I как и раньше идут от 0 до n-1. Поскольку n'' может быть очень большим, массив ''I может переполнить память. Поэтому следует записывать только каждое N''-ое значение аппроксимации. Это возможно с помощью известной нам функции определения остатка: for k in range(1, n+1): ... if k % N 0: # store Итак, каждый раз, когда ''k делится на N'' без остатка, мы записываем значение (в нашем случае каждое сотое). Соответствующая функция представлена ниже. def MCint3(f, a, b, n, N=100): '''Cохраняет каждое N приближение интеграла в массив I и записываем соответствующее значение k' s = 0 I_values = [] k_values = [] for k in range(1, n+1): x = random_number.uniform(a, b) s += f(x) if k % N 0: I = (float(b-a)/k)*s I_values.append(I) k_values.append(k) return k_values, I_values Теперь у нас есть инструмент для того, чтобы посмотреть как изменяется ошибка в интегрировании методом Монте-Карло с ростом n''. Законченная программа выглядит следующим образом, результат работы программы (может несколько отличаться ввиду случайности) представлен справа: thumb|300px|Результат работы программы import random as random_number import matplotlib.pyplot as plt from numpy import array def MCint3(f, a, b, n, N=100): '''Cохраняет каждое N приближение интеграла в массив I и записываем соответствующее значение k' s = 0 I_values = [] k_values = [] for k in range(1, n+1): x = random_number.uniform(a, b) s += f(x) if k % N 0: I = (float(b-a)/k)*s I_values.append(I) k_values.append(k) return k_values, I_values def f1(x): return 2 + 3*x k, I = MCint3(f1, 1, 2, 1000000, N=10000) error = 6.5 - array(I) plt.title('Monte Carlo integration') plt.xlabel('n') plt.ylabel('error') plt.plot(k, error) plt.show() Интегрирование через вычисление площади под кривой Представим, у нас есть некоторая сложная фигура G'' , площадь которой мы хотим вычислить. Эту фигуру мы можем заключить в прямоугольник ''B определяемый координатами X2 x Y2. Далее в этот прямоугольник мы "бросаем" случайные числа из заданного массива Y и считаем сколько точек попало внутрь контура фигуры G''. Далее, зная отношение числа точек, "упавших" в контур ''G к общему числу "брошенных" точек, и умножив это отношение на площадь простой фигуры B'', мы получаем площадь сложной фигуры ''G. Таким же образом мы могли бы находить интеграл от функции, поскольку он и определяет площадь под ней. Представим, мы знаем как выглядит функция, область интегрирования лежит в пределе от a'' до ''b, а максимум функции в указанном интервале равен m''. Тогда аналогично, предыдущему рассмотрению, интеграл будет выражен через число "выброшенных" точек ''N и число упавших "под функцию" M как \frac{M}{N}m(b-a) . Векторизованное решение может в простейшем случае выглядит так: from numpy import random def MCint_area_vec(f, a, b, m, N): x = random.uniform(a, b, N) y = random.uniform(0, m, N) M = y< f(x).size # в квадратных скобках условие, # size считает число ему удовлетворяющих элементов area = M/float(N)*m*(b-a) return area Решение можно написать и через стандартный random в обычном виде с помощью цикла, но скорость его будет существенно меньшей. Категория:Программирование и научные вычисления на языке Python